/*!  \file  MatrixGenerator.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      05/14/2013 11:50:22 AM
 *   \copyright GNU Public License.
 */
#include <math.h>
#include <iostream>
#include <assert.h>
#include <algorithm>
#include <string.h> // memset
#include "Matrix.h"
#include "RandomGenerator.h"
#include "MatrixGenerator.h"
using namespace std;

MatrixGenerator::MatrixGenerator():hasZeroProb(0.23), hasNegativeProb(0.34), hasClustersProb(23),relClusterSigma(1e-5),scale(1.0)  {
	rng = RandomGenerator::defaultGenerator(12343);
}
MatrixGenerator::~MatrixGenerator() {
	if (rng) delete rng; rng = 0;
}

void MatrixGenerator::generate(Matrix & m) const {
	switch (m.getType()) {
		case General        : { generateGeneral    (m); break; };
		case Orthogonal     : { generateOrthogonal (m); break; };
		case Diagonal       : { generateDiagonal   (m); break; };
		case BiDiagonal     : { generateBiDiagonal (m); break; };
		case TriDiagonal    : { generateTridiagonal(m); break; };
		default : { cerr << "MatrixGenerator:: unknown type" << endl;	}
	}
}
void MatrixGenerator::generateGeneral    (Matrix & m) const {
	assert(m.isDense()); // make sure we cna store it
	double * mData = m.getModData();
	if (mData && rng && m.makeGeneral()) 
		rng->fillUniform(m.getStorageSize(), mData);	
		
}
extern "C" {
// LAPACK QR factorization
	void dgeqrf_  ( int * m, int * n, double * A, int * ldA, double * tau, double * work, int * lwork, int * info );
// Generate Q, after QR
	void dorgqr_ ( int * m, int * n, int * k, double * A, int * ldA, double * tau, double * work, int * lwork, int * info );
}
void MatrixGenerator::generateOrthogonal (Matrix & mat) const {
	generateGeneral(mat); 
// perform QR
	int m = mat.rows(), n = mat.cols();
	mat.makeColMajor();	
	double * A = mat.getModData();
	int k        = mat.minDim();
	double * tau = new double[k];\
	double * work = NULL;
	if (tau && A) {
		int lwork = -1;
		int info = 0;
		double wopt = 0;
		dgeqrf_  ( &m, &n, A, &m, tau, &wopt, &lwork, &info );
		lwork = (int) wopt;
		work = new double[lwork];
		if ( work ) {
			dgeqrf_  ( &m, &n, A, &m, tau, work, &lwork, &info );
			if (info == 0) {
				// generate Q after QR	
				dorgqr_  ( &m, &n, &k, A, &m, tau, work, &lwork, &info );
				if (info == 0 )
					mat.makeOrthogonal();
			}
			delete [] work; work = 0;
		}
	}
	if (tau) delete [] tau; tau = 0;
}
void MatrixGenerator::generateDiagonal   (Matrix & m) const {
	if (!m.isSparse()) return;
	bool hasZero     = rng->getUniform() < hasZeroProb;
	bool hasNegative = rng->getUniform() < hasNegativeProb;
	bool hasCluster  = rng->getUniform() < hasClustersProb;
	double min = 0.0, max = 0.0;
	while (max <= min) {
		max = -log( 1 - rng->getUniform());
		min =  exp( rng->getUniform());
	}
	min = hasNegative ? -min : min;
//	cout << endl << "Making matrix " << min << " " << max << endl;
	double * d = m.getModData();
	if (hasZero) d[0] = 0.0;
	if (hasCluster) {
		int numClusters = (int) (m.minDim() * rng->getUniform() * 0.5);
		numClusters = numClusters >= 1 ? numClusters : 1;
		double mean[numClusters];
		for (int i = 0; i < numClusters; i++)
			mean[i] = rng->getUniform(min, max);
		for (int i = hasZero; i < m.minDim(); i++) {
			int c = (int) floor(rng->getUniform() * numClusters);
			assert(c >= 0 && c < numClusters);
			d[i] = rng->getGaussian(mean[c], relClusterSigma * mean[c]); 
		}
		
	} else {
		for (int i = hasZero; i < m.minDim(); i++)
			d[i] = rng->getUniform(min, max); 
	}
	sort(d, d + m.minDim());
}
void MatrixGenerator::generateBiDiagonal(Matrix & bidiag, const Matrix & diag){
	Matrix & m = bidiag;
	double * d = m.getModData();
	unsigned int n = m.minDim();
	assert(bidiag.sameDim(diag));
	memcpy(d, diag.getData(), n * sizeof(double));
	double * e = d + n;
	memset(e, 0, sizeof(*e) * (n - 1));
	for (unsigned int i = 1; i < n; i++) {
		double c,s;	rng->getAzimuth(c,s); // select random cos & sin
		double tmp = 0;
		rot(d[i-1],  tmp, c, s);// apply givens rot from left
		rot(e[i-1], d[i], c, s);
		for (unsigned int j = i; j > 1; j--) {
			getCS(d[j],    tmp, c, s);
			rot(d[j  ],    tmp, c, s); 
			rot(e[j-1], d[j-1], c, s);
			tmp = 0;
			rot(tmp,    e[j-2], c, s);
			getCS(e[j-1],  tmp, c, s);
			rot(e[j-1],    tmp, c, s);
			rot(d[j-1], e[j-2], c, s);
			tmp = 0;
			rot(tmp   , d[j-2], c, s);	
		}	
		getCS(d[1], tmp , c, s); // killEntry(1,0,1, false);
		rot(d[1]  , tmp , c, s);
		rot(e[0]  , d[0], c, s);
	}
}

void MatrixGenerator::generateBiDiagonal (Matrix & m) const {
	if (m.isBiDiagonal() && m.isSquare()) {
		generateDiagonal(m); // fill the diagonal with desired distribution
		double * d = m.getModData();
		unsigned int n = m.minDim();
//		for (int i = 1; i <= n; i++) d[i-1] = i; 
		double * e = d + n;
		memset(e, 0, sizeof(*e) * (n - 1));
//		int step = 0;
		for (unsigned int i = 1; i < n; i++) {
			double c,s;	rng->getAzimuth(c,s); // select random cos & sin
			double tmp = 0;
			rot(d[i-1],  tmp, c, s);// apply givens rot from left
			rot(e[i-1], d[i], c, s);
			for (unsigned int j = i; j > 1; j--) {
				getCS(d[j],    tmp, c, s);
				rot(d[j  ],    tmp, c, s); 
				rot(e[j-1], d[j-1], c, s);
				tmp = 0;
				rot(tmp,    e[j-2], c, s);

				getCS(e[j-1],  tmp, c, s);
				rot(e[j-1],    tmp, c, s);
				rot(d[j-1], e[j-2], c, s);
				tmp = 0;
				rot(tmp   , d[j-2], c, s);	
	
			}	
			getCS(d[1], tmp , c, s); // killEntry(1,0,1, false);
			rot(d[1]  , tmp , c, s);
			rot(e[0]  , d[0], c, s);
		}
	}
//	cout << m;
}
void MatrixGenerator::generateTridiagonal(Matrix & m) const {
	MISSING
}
#ifndef NDEBUG
#include <iomanip>
#include <string.h> // memcpy
#include <math.h>
#include "UnitTester.h"

UNITTEST(MatrixGeneratorTest);
MatrixGeneratorTest::MatrixGeneratorTest() {
	setName 				("Test Matrix Genrator");
	setDescription  	("Test if the generator generats correct matrices");
	setTestAuthor		("Yordan Madzhunkov");
	setDeveloperAuthor("Yordan Madzhunkov");
	setCompileTime    (__DATE__, __TIME__); // each class that inherits this one should also have this line.
};
bool MatrixGeneratorTest::test() const {
	Matrix A(6,6,Orthogonal);
	MatrixGenerator generator;
	generator.generate(A);
	double tolerance = 1e-15;
	bool res = true;
	double norm = (A * A.T()).OrthNorm();
	res &=  norm < tolerance;
	Matrix B(16,10,Orthogonal);
	generator.generate(B);
	norm = (B.T() * B).OrthNorm();
	res &=  norm < tolerance;
	Matrix D(16,16,Diagonal);

	Matrix Bi(6,6, BiDiagonal);
	generator.generate(Bi);
// this code is prototype for generating bidiagonal from diagonal
	double ds[] = {1,2,3,4,5,6};
	Matrix bi(6,6, General);
	Matrix Ds(6,6, Diagonal, ds);
	Matrix Bid(6,6,BiDiagonal);
	generator.generateBiDiagonal(Bid,Ds);

	memset(bi.getModData(), 0, bi.getStorageSize() * sizeof(double));
	bi += Ds;
	for (int i = 1; i < bi.rows(); i++) {
		bi.applyGivensRotationFromLeft(i - 1, i, cos(1.3), sin(1.3));
		for (unsigned int j = i; j > 1; j--) {
			bi.killEntry(j  ,j-1,   j, false);
			bi.killEntry(j-2,j  , j-1, true);
		}
 		bi.killEntry(1,0,1, false);
	}
	return res;
}
MatrixGeneratorTest matrixGeneratorTest;
#endif
